﻿using System;

namespace SoftwareConsulting.BI11.NumericalMethods.Statistics
{
    /// <summary>
    /// Класс хранения статистических моментов случайной переменной
    /// </summary>
    public class StatisticalMoments
    {        
        private double[] _moments; //хранимые моменты

        /// <summary>
        /// Возвращает/устанавливает хранимые моменты
        /// </summary>
        protected double[] Moments
        {
            get
            {
                return _moments;
            }
            set
            {
                _moments = value;
            }
        }

        /// <summary>
        /// Конструктор класса по умолчанию
        /// </summary>
        public StatisticalMoments() : 
            this( 5)
        {
        }

        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="momentsCount">количество хранимых моментов</param>
        public StatisticalMoments(int momentsCount)
        {
            _moments = new double[momentsCount];
            Reset();
        }

        /// <summary>
        /// Сохраняет случайную переменную
        /// </summary>
        /// <param name="x"></param>
        public virtual void Accumulate(double x)
        {
            int n = (int)_moments[0];
            int n1 = n + 1;
            double delta = (_moments[1] - x) / n1;
            double[] sums = new double[_moments.Length];
            sums[0] = _moments[0];
            _moments[0] = n1;
            sums[1] = _moments[1];
            _moments[1] -= delta;
            int[] pascal = new int[_moments.Length];
            pascal[0] = 1;
            pascal[1] = 1;
            double r1 = (double)n / n1;
            int nk = -1;
            n = -n;
            double cterm = delta;
            for (int k = 2; k < _moments.Length; k++)
            {
                sums[k] = _moments[k];
                nk *= n;
                cterm *= delta;
                double term = (1 + nk) * cterm;
                for (int l = k; l >= 2; l--)
                {
                    pascal[l] += pascal[l - 1];
                    term += pascal[l] * sums[l];
                    sums[l] *= delta;
                }
                pascal[1] += pascal[0];
                _moments[k] = term * r1;
            }
        }

        /// <summary>
        /// Возвращает среднее
        /// </summary>
        public virtual double Average
        {
            get 
            { 
                return _moments[1]; 
            }
        }

        /// <summary>
        /// Возвращает количество хранимых моментов
        /// </summary>
        public int Count
        {
            get 
            { 
                return (int)_moments[0]; 
            }
        }

        /// <summary>
        /// Возвращает ошибку в среднем
        /// </summary>
        public double ErrorOnAverage
        {
            get 
            { 
                return Math.Sqrt(Variance / _moments[0]); 
            }
        }
        
        /// <summary>
        /// Возвращает уровень доверия F-теста для указанных статистических моментов
        /// </summary>
        /// <param name="moments"></param>
        /// <returns></returns>
        public double GetFConfidenceLevel(StatisticalMoments moments)
        {
            FisherSnedecorDistribution fisherDistr = new FisherSnedecorDistribution((int)Count, (int)moments.Count);
            return fisherDistr.GetConfidenceLevel(Variance / moments.Variance);
        }
        
        /// <summary>
        /// Возвращает эксцесс (резкость распределения в области максимума)
        /// </summary>
        /// <remarks>Эксцесс нормального распределения, по определению, это ноль</remarks>
        public virtual double Kurtosis
        {
            get
            {
                if (_moments[0] < 4)
                    return double.NaN;
                double kFact = (_moments[0] - 2) * (_moments[0] - 3);
                double n1 = _moments[0] - 1;
                double v = Variance;
                return (_moments[4] * _moments[0] * _moments[0] * (_moments[0] + 1)
                                / (v * v * n1) - n1 * n1 * 3) / kFact;
            }
        }

        /// <summary>
        /// Сбрасывает все моменты
        /// </summary>
        public void Reset()
        {
            Array.Clear(_moments, 0, _moments.Length);
        }

        /// <summary>
        /// Возвращает коэффициент уклона
        /// </summary>
        public virtual double Skewness
        {
            get
            {
                if (_moments[0] < 3)
                    return double.NaN;
                double v = Variance;
                return _moments[3] * _moments[0] * _moments[0]
                                / (Math.Sqrt(v) * v * (_moments[0] - 1)
                                                            * (_moments[0] - 2));
            }
        }

        /// <summary>
        /// Возвращает стандартное отклонение дисперсии
        /// </summary>
        public double StandardDeviation
        {
            get 
            { 
                return Math.Sqrt(Variance); 
            }
        }
        
        /// <summary>
        /// Возвращает уровень доверия T-теста для указанных статистических моментов (дисперсия множеств может отличаться, для этого и нужно приближение)
        /// </summary>
        /// <param name="moments"></param>
        /// <returns></returns>        
        public double GetTApproximateConfidenceLevel(StatisticalMoments moments)
        {
            StudentDistribution studentDistr = new StudentDistribution(
                                                (int)(Count + moments.Count - 2));
            return studentDistr.GetConfidenceLevel((Average / StandardDeviation
                                                - moments.Average
                                                / moments.StandardDeviation)
                                                / Math.Sqrt(1 / Count
                                                            + 1 / moments.Count));
        }

        /// <summary>
        /// Возвращает уровень доверия T-теста для указанных статистических моментов (дисперсия обоих множеств должна быть одинаковой)
        /// </summary>
        /// <param name="moments"></param>
        /// <returns></returns>
        public double GetTConfidenceLevel(StatisticalMoments moments)
        {
            int dof = (int)(Count + moments.Count - 2);
            double sbar = Math.Sqrt((UnnormalizedVariance + moments.UnnormalizedVariance) / dof);
            StudentDistribution studentDistr = new StudentDistribution(dof);
            return studentDistr.GetConfidenceLevel((Average - moments.Average)
                                            / (sbar * Math.Sqrt(1.0 / Count
                                                            + 1.0 / moments.Count)));
        }

        /// <summary>
        /// Возвращает ненормализированную дисперсию
        /// </summary>
        /// <remarks>необходимо для T-теста</remarks>
        public virtual double UnnormalizedVariance
        {
            get 
            { 
                return _moments[2] * _moments[0]; 
            }
        }

        /// <summary>
        /// Возвращает дисперсию
        /// </summary>
        /// <remarks>дисперсия включает в себя коррекционный фактор Бесселя</remarks>
        public virtual double Variance
        {
            get
            {
                return _moments[0] < 2
                        ? double.NaN
                        : UnnormalizedVariance / (_moments[0] - 1);
            }
        }
    }
}
